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ordinary differential equations over continuous variables representing 
the concentrations of the reactant species. It is known that two net- 
works may possess the same governing mass-action dynamics despite 
disparate network structure. To date, however, there has only been 
limited work exploiting this phenomenon even for the cases where one 
network possesses known dynamics while the other does not. In this 
paper, we bring these known results into a broader unified theory which 
we call conjugate chemical reaction network theory. We present a the- 
orem which gives conditions under which two networks with different 
governing mass-action dynamics may exhibit the same qualitative dy- 
namics and use it to extend the scope of the well-known theory of 
weakly reversible systems. 



1 Introduction 

A chemical reaction network is given by sets of reactant species reacting at 
prescribed rates to form sets of product species. Under appropriate assump- 
tions (well-mixing, constant external conditions, large number of reacting 
molecules, mass-action kinetics, etc.) such systems can be modeled deter- 
ministically by an autonomous set of ordinary differential equations over 
continuous variables representing the concentrations of the reactant species. 
The resulting mathematical models have a long history and enjoy applica- 
tions in fields such as systems biology, industrial chemistry, pharmaceutics, 
etc. [5| ll5ll22"| l25j. Active theoretical topics within the mathematical litera- 
ture include determining network conditions which guarantee locally stable 
dynamics [1401801902611270321 . conditions which prohibit/permit multistabil- 
ity 0E1CLIJG1] > an d conditions which ensure the persistence of the chemical 
species [HElElElEn] . 

One topic which has gained only fleeting attention until recently has 
been determining conditions under which networks with different reaction 
structure exhibit the same qualitative dynamics. The most comprehensive of 
these results is Craciun and Pantea [10] wherein the authors consider condi- 
tions under which a mass-action system does not permit uniquely identifiable 
rate constants in the sense that two disparate reaction networks could pro- 
duce the same governing set of differential equations under the assumption 
of mass-action kinetics for specified choices of the rate constants. Within 
the literature, this is sometimes referred to as the "fundamental dogma of 
chemical kinetics" |12 p29ll30] . This work has recently been continued by G. 
Szederkenyi et al. who have presented computational methods for deter- 
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mining networks which have equivalent dynamics [29H31] . Related work in 
the literature can be found in [^ [T2l[T9H2"T] . 

In this paper, we attempt to bring these results together into a unified 
framework and language. We have chosen to borrow from dynamical systems 
theory in calling two mass-action systems with related dynamics conjugate 
systems [23, 33 . More specifically, we will call two mass-action systems 
conjugate if there is a mapping which takes trajectories of one system into 
trajectories of the other. We will focus on the subset of linearly conjugate 
systems. To the best of our knowledge, all of the results contained in the 
literature to date on qualitative equivalence of mass-action systems restrict 
their focus to conditions under which the governing differential equations for 
two networks are identical; consequently, the required conjugacy mapping 
is the identity. We will go further than these results with Theorem 13.21 in 
giving conditions on the network for which a non-trivial linear mapping is 
required to demonstrate conjugacy. 

In practice, we are typically presented with a specified network and asked 
to determine the dynamic behaviour. Linear conjugate chemical reaction 
network theory presents a powerful approach to solving such a problem: if 
a network with unknown dynamics can be shown to be linearly conjugate 
to a network with known dynamics, the dynamical properties of the second 
apply to the first as well. While the second network is generally unknown, 
the apparent dynamics of the first network may provide clues as to the 
class of networks within which to search for a linearly conjugate network. 
We focus particularly on the class of weakly reversible networks considered 
in [14 y i8 1 fT9] . Szederkenyi and Hangos have considered some aspects of this 
problem from a computational perspective |30j . 

The paper is organized as follows: in Section [2] we introduce the termi- 
nology and definitions which will be used throughout the paper; in Section [3] 
we introduce our notion of conjugate chemical reaction networks, present the 
relevant results from the literature, give a new result (Theorem 13. 2p which 
goes beyond this known theory, and connect our results to weakly reversible 
networks through examples; we close in Section [J] with some reflections on 
the results obtained and avenues for future work. 

Throughout the paper, we will let M> and M> denote the m-dimensional 
spaces with all coordinates strictly positive and non-negative, respectively. 
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2 Background 



In this section, we introduce the concepts and notation relevant to the study 
of chemical reaction networks. We also summarize the work conducted to 
date on the subject of conjugacy of chemical reaction networks. 

2.1 Chemical Reaction Networks 

We will let Aj denote the species or reactants of the network and define 
\S\ = m where S is the set of distinct species of the network. A chemical 
reaction network consists of the union of the elementary reaction set 

m m 

IZi : ^ z ij^j y ^ Z ij-Aj, i = 1, . . . , r (1) 

3=1 3=1 

where Zij,z[j G Z>o are the stoichiometric coefficients of the i th reaction. 
The set of all reactions in the network will be denoted 1Z so that \H\ = r. 
It will occasionally be convenient to adopt the short-hand Cj = YHj=\ z ij^j 
and C[ = Y^IjLx Z ij-Aj so that the reaction network can be given by 

TZi-. d — > C'i i = l,...,r. (2) 

In this setting, d will be referred to as the reactant complex and C[ will 
be referred to as the product complex. A chemical reaction network will be 
denoted by the triplet M = (S,C,1Z). 

There is no demand that the stoichiometry of each reactant and product 
complex be distinct except that Cj ^ C\ for i = 1, . . . , r (so each reaction is 
meaningful) . That is to say, a single stoichiometrically distinct complex may 
be the reactant or product complex for numerous reactions. It will occa- 
sionally be convenient to index the complexes not according to the reactions 
they occur in, but by their stoichiometric distinctiveness. In this setting, 
reactions are not represented as a list but as interactions between the n 
stoichiometrically distinct complexes of the network. We will denote by C l , 
i = 1, . . . , n, the n stoichiometrically distinct complexes of the network and 
let C denote the set of such complexes. 

We are often interested in the algebraic structure of a reaction network. 
These properties will be most readily seen by indexing the complexes by their 
stoichiometric distinctiveness according to C l , i = l,...,n. The following 
concepts will be of interest to us. 
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Definition 2.1. Consider a chemical reaction network N . We will say there 
is a path from C l G C toC 3 G C if there exists a sequence of stoichiometrically 
distinct complexes {C Ul , C U2 , . . . , C u <} , such that C* = C 1 , C 3 = C Ul , and 

C Vk — > C Uh+1 G n 

for every k = 1, . . . ,1 — 1 . 

Definition 2.2. A chemical reaction network M is said to be weakly re- 
versible if the existence of a path from C l to C 3 implies the existence of a 
path from C 3 to C l . 

Intuitively, weak reversibility means every reaction can be reversed by 
a suitable sequence of subsequent reactions. It is not necessary that every 
reaction have a reverse reaction (i.e. that the system be reversible), as can 
be seen by the weakly reversible network 

C 1 — ► C 2 

\ y 

c 3 . 

2.2 Mass- Action Kinetics 

We are particularly interested in the evolution of the concentrations of the 
chemical species as time progresses. We will let x\ = [At] denote the con- 
centration of the i th species and denote by x = [x± X2 • • • x rn ] T G M> the 
concentration vector. 

The differential equations governing the chemical reaction network J\f 
given by (jTJ under the assumption of mass-action kinetics are given by 

^ = f (x) = J> ft - *) X" 1 (3) 

i=l 

where x Zi = IIj=i x /* an( ^ ^« > is the rate constant for the i th reaction. 
We will denote by k = [k\ &2 • • • kr] T G M> ^he rate constant vector. 

Chemical reaction networks endowed with mass-action kinetics (j3J) will 
be called mass-action systems and will be denoted by (S, C,1Z, k). Since 
mass-action kinetics is the only form of kinetics considered in this paper, we 
will use M interchangeably to denote both the reaction network (S,C,H) 
and the mass-action system (S,C,H, k). 

Several fundamental properties of mass-action systems are readily seen 
from (j3J). In particular, it is clear that solutions are not able to wander 
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around freely in K™ . The following concepts help clarify where solutions 
may lie. 

Definition 2.3. The kinetic subspace for a mass-action system is the 
smallest linear subspace S* C M m such that 

S* = im(f(x)) . 

The dimension of the kinetic subspace will be denoted by \S*\ = s* . 

Definition 2.4. The stoichiometric subspace for a mass-action system 
is the linear subspace S C R m such that 

S = span { — Zj) | i = 1, . . . , r} . 

The dimension of the stoichiometric subspace will be denoted by \S\ = s. 

While the stoichiometric and kinetic subspaces are related, they need not 
coincide, in which case we have the strict inclusion S* C S (see Example 2 
in Section I3.3P . Since we clearly have f (x) G S and f (x) G S* , by integrating 
([3]) we can see that solutions can be partitioned according to their initial 
conditions xo G K> m ^ ne following way. 

Proposition 2.1 ( [T91I32] ). Let x(t) be the solution to 0j with x(0) = x G 
K^ . Then x(t) G C* = (x + 5*) n for t > w/iere C* is called the 
positive kinetic compatibility class associated with xq G 

In other words, the sets C* Q are invariant spaces of the dynamics ©. 
A more commonly encountered invariant space of (|3|) is the positive stoi- 
chiometic compatibility class C Xo = (xo + S) Pi K> - ^he sets C Xo and C* Q 
are related by the following result, which allows us to restrict attention to 
C Xo for weakly reversible networks. 

Lemma 2.1 (Corollary 1, [E]). For weakly reversible mass-action systems 
the stoichiometric and kinetic subspaces coincide. 

2.3 Complex Balanced Systems 

Many mass-action systems can be classified according to the nature of the 
equilibrium concentrations they permit. A particularly important class of 
such systems, which will be considered extensively in this paper, are complex 
balanced systems (see ) ■ 
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Definition 2.5. An equilibrium concentration x* £ ]R> 0/ a mass-action 
system ^ will be called complex balanced if, for every stoichiometrically 
distinct complex C° E C, we have 

£ fc,(x7* = (x«)- £ 

j=i i=i 
c<=c° C J= C° 

^4 mass-action system will be called a complex balanced system if every 
equilibrium concentration is a complex balanced equilibrium concentration. 

In other words, an equilibrium concentration is complex balanced if the 
net flow rate into every stoichiometrically distinct complex is exactly bal- 
anced by the net flow rate out of the complex. It is known that every 
mass-action system permitting a complex balanced equilibrium concentra- 
tion corresponds to a weakly reversible reaction network (Theorem 2B, [18]). 

Analysis of complex balanced systems is made easier by the following 
lemma. 

Lemma 2.2 (Lemma 5B, |19j). If a mass- action system is complex bal- 
anced at some concentration x* 6 ^>o> then it is complex balanced at all 
equilibrium concentrations. 

Consequently, any mass-action system permitting a complex balanced equi- 
librium concentration is a complex balanced system. It is worth noticing, 
however, that a mass-action system may be complex balanced for some val- 
ues of the rate constants and not for others. 

2.4 Locally Stable Dynamics 

Mass-action systems permit a variety of behaviours, including multistability, 
periodic behaviour, and chaotic behaviour (see [13]). In this paper, we will 
be particularly interested in networks exhibiting locally stable dynamics. 

Definition 2.6. A mass-action system is said to exhibit locally stable 
dynamics if there is a unique positive equilibrium concentration x* € M™ 
within each kinetic compatibility class C* and that equilibrium concentration 
is locally asymptotically stable relative to C* . 

Locally stable dynamics is a particularly desirable form of dynamics since 
the behaviour is so regular and predictable; however, it is local in nature and 
does not preclude trajectories from approaching <9C* Q . Mass-action systems 
for which each x* is a global attractor for C* are called globally stable. The 
following result was derived by Horn and Jackson in [19], 
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Theorem 2.1 (Theorem 6A and Lemma 4C, [19J). If a mass-action system 
is complex balanced then it possesses locally stable dynamics. 

Since complex balanced systems are necessarily weakly reversible, it follows 
by Lemma 12.11 that we may restrict our attention to the positive stoichio- 
metric compatibility classes C X[) in the definition of locally stable dynamics 
for such systems. 

This result is powerful in that it guarantees a particularly desirable form 
of dynamics based solely on consideration of the equilibrium set of the sys- 
tem. Together with Feinberg, the authors go further in relating the complex 
balancing condition to the algebraic structure of the reaction network ([2]) . In 
particular, they demonstrate that every weakly reversible mass-action sys- 
tem is complex balanced for at least some choice of the rate constants [14018] . 
They also show that there are algebraic conditions on a weakly reversible 
network ([2]) sufficient to guarantee the associated mass-action system is com- 
plex balanced for all rate constants (Theorem 4A, |18|). Further consider- 
ation of what the conditions on the rate constants required to guarantee 
complex balancing look like is contained in [6]. 

3 Conjugate Chemical Reaction Networks 

The focus of this paper is determining when two mass-action systems can be 
shown to have qualitatively identical behaviour despite disparate reaction 
networks (pQ). Our approach to this problem is to show that there is a 
suitably nice mapping between the flows of ([3]) for one network and another. 
In the standard theory of differential equations, this notion is captured in the 
well-studied concept of conjugacy (for a complete introduction to conjugacy 
of dynamical systems, see [231133] ) . 

We define the reaction network M' = (S,C',1Z') as consisting of the 
reaction set 

m 

3=1 

or, alternatively, 

Tl'i : Ci — ► C[, i = l,...,f. (5) 

The associated mass-action system will be denoted (S,C ,TZ' , k') where the 
rate constants are given by k\ > 0, i = 1, . . . , r. We will let <l?(xo,t) denote 
the flow associated with the mass-action kinetics ([3]) for N and ^(yo,^) 



1 f 

X, . . . , I 



(4) 



3=1 
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denote the flow associated with the mass- action kinetics © for JV 7 . We will 
adopt the convention of referring to M as the original network and N' as 
the target network. 

Note that, while we follow the notation of |10] in denoting any second 
network by A/ 7 = (S,C',1Z'), we distinguish the relevant components of the 
second network using tildes to avoid confusion with the vectors from the 
first system. Also notice that the networks W and A/"' are allowed to have not 
only different complexes and reactions, but different numbers of complexes 
and reactions; the number of species \S\ = m, however, is required to be the 
same. 

We are now prepared to define our notion of conjugacy of chemical re- 
action networks. 

Definition 3.1. Consider two mass-action systems N and N' . We will say 
M andM' are C k -conjugate if there exists a C k -diffeomorphism h : M> i— 
M™ such that h($(x ,i)) = *(h(x ),i) for all x G R™ . 

Definition 3.2. We will say M and N' are linearly conjugate if they are 
C°° -conjugate and the diffeomorphism h : M™ i— > ]R> is linear. 

In this paper, we will focus on the notion of linear conjugacy. Note that 
any linear diffeomorphism is necessarily C°° so that any linear conjugacy is 
a C°°-conjugacy. The following results clarify the form the linear mapping 
h : 1R™q I—?- IR^q may take and the implications of conjugacy. 

Lemma 3.1. A linear, bijective mapping h : R™ i— > R™ may consist of at 
most positively scaling and reindexing of coordinates. 

Proof. Consider a linear, bijective mapping h : M> i— >■ M> . Since h(x) is 
linear, it can be represented h(x) = Ax where A G jj mxm an( j s j nce h(x) 
is bijective, it has an inverse h _1 (x) = A~ l x. Since the mappings are from 
M™q to K>o> au entries m A and A -1 must be non-negative and every row 
of A and A" 1 must contain at least one non-zero entry. 

Suppose there is a row of A with more than one non-zero entry. Since A 
and A~ l may contain no negative numbers, in order to satisfy A A~ l = I 
this implies that there are at least two rows of A -1 which contain zeroes 
in the same m — 1 columns. Such an A , however, would have a zero 
determinant and therefore be non-invertible, which is a contradiction. 

It follows that each row of A has precisely one positive entry. Since A 
is invertible it follows that each column of A also has precisely one positive 
entry so that A is a positively weighted permutation matrix. In terms of 
the transformation h(x) = Ax this means the mapping may only positively 
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scale and re-index the components of the vector x, which completes the 
proof. □ 



Lemma 3.2. If a mass-action system N is linearly conjugate to a mass- 
action system N' and Af' exhibits locally stable dynamics, then N exhibits 
locally stable dynamics. 

Proof. The result follows trivially from Lemma 13. II and Definition 13.11 □ 

It is worth noting that other qualitative properties of mass-action sys- 
tems are also preserved by linear conjugacy (multistability, persistence, 
boundedness, etc.). Some aspects of qualitative equivalence of M and M 1 
can, however, fail for non-linear conjugacies (see the discussion in Section U 
and the Appendix). 

3.1 Known Results 

In this section, we give a brief summary of the results which are, to the 
best of our knowledge, the only attempts to demonstrate conjugacy of two 
mass-action systems. In all of the results described in this section, con- 
jugacy is demonstrated by showing an exact equivalence between the gov- 
erning differential equations ([3]) for M and M 1 . (This phenomenon is called 
macro- equivalence in [19] and confoundability in [10] , In [29] and the related 
literature, two networks with identical dynamics are called two realizations 
of the same reaction kinetic differential equations.) 

The most thorough study of conjugacy to date has been conducted by 
Craciun and Pantea [10]. In that paper, the authors considered the problem 
of experimentally assigning values to rate constants to systems with lin- 
early dependent reactions flowing from the same reactant complexes. They 
present the following result. 

Theorem 3.1 (Theorem 4.4, [10]). There exist rate constants choices such 
that the mass- action systems M and N' are conjugate with h(x) = x if and 
only if they have the same reactant complexes and Cn(C°)PiC'ji>(C ) ^ for 
every reactant complex C°, where 

C-ji(C°) = l^^a^z^ — Zj) | Qj > if Zi = z°,Qj = oifterwisej' . (6) 

This theorem gives necessary and sufficient conditions for two chemical 
reaction networks M and N' to admit rate constant vectors which generate 
the same set of governing differential equations ([3]). It is clear that conjugacy 
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follows according to Definition 13.11 It should be noted, however, that two 
reaction networks may be conjugate for some choices of the rate constants 
and not for others. (It was noted in [28] that this result is deficient in 
that it ignores the possibility that the net flow from a reactant complex 
in either M or A/ 7 could equal zero. Consequently, conjugacy could hold 
for networks with different reactant complexes so long as the corresponding 
outflows cancel in Q. We will ignore this sensitivity for the time-being.) 

In practice, we are typically asked to determine the dynamics of a mass- 
action system N without reference to any potential target networks N' . 
The results of [10] propose no methodology to find a potential JV 7 satisfying 
Theorem 13. 1[ i.e. the target network has to be specified. Szederkenyi has 
extended the work contained in |10j by presenting computer algorithms ca- 
pable of finding networks N' with equivalent mass-action kinetics ([3]) as N . 
He calls the networks N and N' alternative realizations of the kinetics ([3]) . 
In [29] he gives an algorithm for finding sparse and dense realizations (i.e. 
realizations with the fewest and greatest number of reactions) and in [30] 
he and Hangos give an algorithm for determining detailed and complex bal- 
anced realizations for specified kinetics ([3]). 

Averbukh also considers conditions which relate the dynamics of an un- 
determined network J\f to a network N' with known dynamics. In particular, 
he presents conditions under which a general network has the same dynam- 
ics as a detailed balanced network (Theorem 2 of [4]). This is a powerful 
result since detailed balanced networks are known to exhibit locally stable 
dynamics (151152] . 

Conjugacy is also considered by Krambeck in Section 6 of |20] for de- 
tailed balanced systems where it is referred to as non-uniqueness of the rate 
constants. Horn and Jackson briefly consider conjugate systems in their 
seminal paper |19j . Their primary example is the network 

2Ai + A 2 3Ai 

3A 2 A- A! + 2A 2 

where e > 0. They show that the network exhibits locally stable dynamics 
for e > 1/6 and that the network possesses the same mass-action kinetics 
([3]) as a complex balanced network M' for e > 1/2. The Master's thesis of 
MacLean also contains specific examples of networks which are conjugate 
to complex balanced systems |21j . This connection with complex balanced 
systems is made more explicit in her unpublished research notes. 
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3.2 Original Results 

In this section, we present our main original result regarding conjugacy 
of chemical reaction networks. In Section 13.31 we show how Theorem 13.21 
can broaden the scope of weakly reversible networks theory through several 
illustrative examples. 

In what follows, we will let C Teac t denote the set of reactant complexes 
in either the complex set C or the complex set C 

Theorem 3.2. Consider two mass-action systems N and N' . Suppose that 
for the rate constants ki > 0, i = l,...,r, there exist constants hi > 0, 
i = 1, . . . , f, and Cj > 0, j = 1, . . . , m, such that, for every C° € C reac t, 



ki«-*i)=T UZ'i-ii) (7) 

i=l i=l 
Ci=C° Ci=C° 

where T =diag{cj}J l =l . Then M is linearly conjugate to M' with rate con- 
stants 

in 

h = bij{cf, i = l,...f. (8) 

j=i 

It is important to note that the reactant complex set for C need not be 
the same as that of C. When C° G C reac t is not an element of the reactant 
complex set of C, we will consider the summation on the left-hand side of 
([7]) to be empty, and similarly for the right-hand side of ([7]) when C° is not 
an element of the reactant complex set of C, i.e. 

r f 

Y hi^'i - Zj) = and ^ ^(z- - z;) = 0, 

i=l i=l 
C t =C° c,=C0 

respectively. In order to satisfy ([7|), therefore, if one system contains a 
reactant complex not contained in the other, it is necessary for the origin to 
lie in the cone generated by the reaction vectors flowing from that reactant 
complex in the other system. 

Proof. Let ^(xo, t) correspond to the flow of the mass-action system ([3]) asso- 
ciated to the reaction network N given by (pQ) . Consider the linear mapping 
h(x) = T _1 x where T =diag{c ? -}™ =1 . Now define ^(yo,t) = T~ 1 <I > (xo,t) so 
that $(x ,t) = T*(y ,t). 
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Since <3?(xo,i) is a solution of ([3]) for the reaction set (P), we have 



tt'(y ,t) = T- 1 cJ»'(x ,t) 



r 



= T- 1 ^^(z^-z i )$(x ,t) Zl 



i=l 



r 



= T_1 E E fci«-«i)*(xb.*)"° 



c°ec reQ rf »=i 



C i= C° 



r 



Ci=C° 



r m 



It is clear that Vl/(yo,i) is the flow of ([3]) for the reaction network with 
rate constants given by ([5]). We have that h($(xo,t)) = \£(h(xo), i) for all 
x £ M> and t > where y = h(x ) since y = ^(yo, 0) = T~ 1 $(x , 0) = 
T _1 xq. It follows that the networks J\f and J\f' are linearly conjugate by 



This result gives conditions under which two mass-action systems N 
and M' are linearly conjugate. This is a particularly useful result when the 
qualitative properties of the original network are obscure while the behaviour 
of the target network is well understood. 

With the exception of the scaling matrix T, condition ([7]) is very similar 
to the cone intersection condition in Theorem 13 . 1 1 where the constants hi > 0, 
i = 1, . . . , r, and bi > 0, i = 1, . . . , f, correspond to the magnitudes of the 
cone generators (i.e. the reaction vectors). If we allow ki and bi to vary we 
have 



Definition 13.21 and we are done. 



□ 




and 
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according to ©. 

The following two results can be obtained from Theorem 13.21 by al- 
lowing the rate constant vector k £ MJ >0 to vary. In these results we let 
T =diag{cj}™ =1 and consider Cj > 0, j = 1, . . . , m, to be fixed. 

Corollary 3.1. Consider two mass-action systems M and N' . Then there 
exist rate constant vectors k £ M> and k' £ M> suc/i i/iai M and M' are 
linearly conjugate with h(x) = T _1 x if and only if for every C° £ C reac t we 
have C n {C°) n [T C^(C )] + 0. 

Corollary 3.2. Consider two mass- action systems N and A/ 7 . T/ieri /or 
every rate constant vector k £ M> there exists a rate constant vector k' £ 
R> suc/i that M is linearly conjugate to M' with h(x) = T _1 x if and only 
if for every C° £ C react we have Cn{C°) C [T C^/(C )]. 

Proof. The forward implications follow directly from Theorem 13.21 

To prove the 'only if portions of the results, notice that the assumption 
of conjugacy with h(x) = T _1 x implies that 

r I m \ 

*W) = E E T_1 II c ? « - z *) *(yo. *) zo (9) 

C°eC re[lc t i=l \ i=i / 

C l= C° 

while we have 

f 

*W) = E E ^-^(yo,*) 20 (10) 

C GC reac t i=l 
C l= C° 

from ([3]). In order to have equality between ([9]) and (jlpp we require that 

f 

i=i 
C l= C° 

for every C° £ C reac t- The desired cone conditions follow immediately from 
the conditions on the rate constants vectors k £ RC_ and k' £ W >Q . □ 

3.3 Weakly Reversible Networks 

In Section [3l the results depended on having two given networks M and A/ 7 
to compare. In standard practice, however, we have only a single network 



=1 
=e<> 
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M whose dynamics are unknown and we need to find the target network J\f' 
whose dynamics are understood. 

In this section, we will consider a particularly broad and well-understood 
class of such target networks in weakly reversible networks. Since it is known 
that weakly reversible systems are complex balanced for at least some values 
of the rate constants, and therefore exhibit locally stable dynamics for those 
rate constants values, it is a reasonable starting point when considering a 
network M which seems to exhibit locally stable dynamics to search for 
weakly reversible target networks M 1 to which it could be conjugate. 

In practice, however, there are many sensitivities which can arise in 
choosing a suitable target network N' which is weakly reversible. We will 
illustrate the applicability, and limitations, of Theorem 13.21 to such cases 
through four examples. The first is an example where linear conjugacy 
to a weakly reversible network which exhibits locally stable dynamics can 
be universally shown. The second is an example where linear conjugacy 
to a weakly reversible network can only be shown for certain choices of 
the rate vector k £ W >Q . This is also an example where S* and S do 
not always coincide for the original network N . The third is an example 
where linear conjugacy to a weakly reversible network holds universally but 
conditions on k S BL> are still required to guarantee locally stable dynamics. 
This example also demonstrates how these conditions can be reduced by 
creatively "splitting" a reaction in the target network A/ 7 . The fourth is an 
example where a "phantom" reactant complex is required to demonstrate 
linear conjugacy to a weakly reversible network. 

Our general technique in this section will be to search for weakly re- 
versible target networks N' which involve the same reactant complexes as 
the original network N . 

Example 1: Consider the chemical reaction network M given by 

M _ Ax + 2A 2 A Ax + ZA 2 Ai + A 2 3Ai 

2Ai A A 2 - 

We want to find a weakly reversible target network N' involving the 
reactant complex set of the original network M, which is 

{Ai + 2A 2 ,A 1 + 3A 2 ,Ai+A 2 ,2Ai}. 

Many such networks can be eliminated for failing to be weakly reversible, 
leaving a relatively small set of possibilities. One such possibility is given 
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by 



TV' : 



Ai + 2A 2 t±Ai + 3A 2 
A1 + A2 ^ 2Ai. 

k<4 



In order for TV and TV' to be conjugate, we need to find bi > 0, Cj > 0, 
z = 1, ... ,4, j 



1,2, such that 



fci 




1 



h 

k 1 



1 



bi 

b 2 
b, 
bi 



a 

c 2 

a 

c 2 

ci 

c 2 

a 

c 2 




1 



1 



k 3 , 64 



fe 4 , Cl 



It can be easily verified that b\ = k\, 62 = 2/c 2 , 63 
c 2 = 1 works. 

It follows from ([8]) that TV is linearly conjugate to TV' with rate constants 
given by k\ = 2k\, /c 2 = 4/c 2 , /C3 = 2k$ and /C4 = 4^. It is known that 
TV' is complex balanced, and therefore exhibits locally stable dynamics, for 
every rate constant vector k' E M>o- It follows that the original network TV 
exhibits locally stable dynamics for every rate constant vector k £ ]R> by 



Lemma [ 

It could also be noted that C n (C°) = T C^/(C°) for every G C reac t so 
that TV and TV' satisfy the hypotheses of Corollary 13.21 and therefore linear 
conjugacy holds unconditionally. 

Example 2: Consider the chemical reaction network TV given by 



TV : 



A 2 



— A! A- 2^ 2 -^3Ai. 



The only weakly reversible target network involving the same reactant 
complex set as TV is 

TV' : Ai S 2^2- 

k 2 
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In order to satisfy (J7J), we need to find b\ > 0, 62 > 0, c\ > 0, and c 2 > 
such that 



fci 

h 

while satisfying (|8|) requires 
fei = hex, 



' -1 " 




1 


= 61 


1 


+ h 


-2 



Cl 







C2 



Cl 







C2 



and 



b 2 (? 2 . 



(12) 



The system (fTTj) corresponds to satisfying &i = 61 ci, fci = 261C2, k 2 + 
3k$ = 62 ci, and &2 + ^3 = ^2C2- From the first two equations, we have 
C1/C2 = 2 while from the last two we have c\/c 2 = (k2+Sks)/(k 2 +ks), which 
implies ([lip can be satisfied if and only if k 2 = k^. With this restriction, 
the system can be satisfied for b\ = k±, 62 = 2k 2 = 2k%, c± = 2, C2 = 1. It 
follows from (|12p that fei = 2k\ and A;2 = 2/^2 = 2k^. 

It is known that M 1 is complex balanced, and therefore exhibits locally 
stable dynamics, for all values of k\ > and k 2 > 0; however, because we 
required a condition on the rate constants of M in order for condition (jlip 
to be satisfied, ftf does not exhibit locally stable dynamics unconditionally. 
In fact, it exhibits locally stable dynamics only for k 2 = k^. For k 2 > ks, 
all trajectories tend to the origin, while for k% > k 2 all trajectories become 
unbounded. 

It is worth noting that the kinetic subspace S* is two-dimensional for 
J\f for all rate constants values except k 2 = k% when it collapses to a single 
dimension and we have the strict inclusion S* C S. Since N' is weakly 
reversible, we always have S* = S for N' by Lemma 12. II and we notice that 
this is always one-dimensional. The systems will only be conjugate when 
the dimensions of the kinetics compatibility classes match, which only occurs 
when k 2 = k%. 

It could also be noted that C n {C )n[T C^iC )} / for every C° G C react 
but Cti{C°) £ [T C n i{C )} for C° = 2A 2 G C react . Consequently, the net- 
works J\f and J\f' satisfy the hypotheses of Corollary 13.11 but not Corollary 
13.21 linear conjugacy with h(x) = T _1 x cannot therefore be guaranteed for 
all rate constant vectors k G K>q- 
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Example 3: Consider the chemical reaction network N given by 



TV 



Ai + 2A 2 
2Ai + A 2 
Ai + 3^2 



A l 
2>A 2 

Ai+A 2 



3Ai + A 2 



for e > 0. 

We search for target networks N' with reactions flowing between the 
complexes in the reactant set of J\f, which is 

{Ai + 2A 2 , 2Ai + A 2 , A! + 3A 2 ,Ai + A 2 } ■ 

Many such networks can be eliminated for failing to be weakly reversible, 
leaving a relatively small set of possibilities. We will choose the network N' 
given by 



Ai + 2A 2 A Ai + A 2 



&4 -j^ &5 

Ai + 3^2 < 



k 2 

2Ai + A 2 



where we have chosen to "split" the reaction flowing from the reactant com- 
plex Ai + 3^2 into two weighted reactions. The utility of this technique will 
become apparent momentarily. 

In order to satisfy ([7]) we need to find constants b{ > 0, Cj > 0, i = 

5, j 



1,2, such that 







-2 

2 


-2 



2 




b 2 



ci 
c 2 

ci 
c 2 

ci 
c 2 

ci 
c 2 







1 





+ b 5 



We will choose the solution set bi = 2e, b 2 



1, h = 2(1 - t), h = t, 



ci = 2, c 2 = 1 where < t < 1 is a weighting constant. This gives rise to the 
rate constants ki = 4e, k 2 = 2, = 4, k^ = 4(1 — t) and k§ = 2t according 
to (ED. 
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There is one condition on the rate constants of M 1 in order for the mass- 
action system to be complex balanced. That condition is 

1 - t 
€ ~ y/2=t' 

Each < t < 1 corresponds to specific network M' which is conjugate to N . 
Consequently, we can guarantee that M is conjugate to a complex balanced 
system, and therefore exhibits locally stable dynamics, for the range of values 
< e < 1/ y/2. Notice that if we had not split the reaction flowing from the 
complex A\ + 3A2 and instead had all of the weight represented in k^, we 
would have only been able to show that M exhibits locally stable dynamics 
for e = 1/V2. 

It can be checked empirically that locally stable dynamics appears to be 
exhibited for N for all values e > 0. There is a unique positive equilibrium 
concentration given by (xl,x* 2 ) = (1, -e/2 + + 1) which is locally 

asymptotically stable for all e > according to standard linearization the- 
ory. It is not our claim, therefore, that this theory represents a complete 
classification of locally stable dynamics, even for networks which exhibit lo- 
cally stable dynamics for some values of the rate constants. 

Example 4: Consider the chemical reaction network J\f given by 

Ai 2A l + 2A 2 A A 2 Ai+A 2 . 

A quick analysis ([3]) reveals that N appears to exhibit locally stable 
dynamics which suggests that the network may be linearly conjugate to a 
complex balanced system. It can be quickly seen, however, that there is no 
weakly reversible network A/ 7 involving the same reactant complexes as M 
which satisfies the requirements of Theorem 13.21 

Theorem 13.21 however, does not require the target network M 1 to use the 
same reactant complexes as M. We can have a reactant complex C° from 
the reactant complex set of C which is not in the reactant complex set of C 
so long as 

f 

E b ^ - ^) = °- 

i=l 
C>=C» 

One possible target network M' which makes use of such a "phantom" re- 
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actant complex is given by 



M' : 



A1+A2 



A 2 



A! 




\ f 2 

2Ax + 2A 2 . 



Since the reactions corresponding to k\, k 2 and k 3 in M are the same as 
those for k±, k 2 , and k^ in J\f, we set c\ = c 2 = 1, b\ = k±, b 2 = k 2 , and 
^3 = k^. Since the complex A\ +A 2 is not a reactant complex in C, in order 
to satisfy ([7D it is required that 



ki 








" 1 " 




' -1 " 




" " 


-1 


+ ^5 


1 


+ h 











This can be satisfied if k^ = k§ = k§ = t for any t > 0. 

There is one condition on the rate constants of N' required for the net- 
work to be complex balanced. That condition in terms of the rate constants 



Since every value of t > corresponds to a valid conjugate network A/ 7 , this 
is no restriction at all. It follows that the original network M is linearly 
conjugate to a complex balanced system for all choices of rate constants and 
therefore universally exhibits locally stable dynamics by Lemma 13.21 

In other words, we are able to demonstrate the network M is linearly 
conjugate to a complex balanced network M' , and therefore possesses lo- 
cally stable dynamics, by adding a "phantom" reactant complex which con- 
tributes no dynamical information to the mass-action kinetics ([3]) . It should 
be noted, however, that in order for the target network N' satisfying (|7|) to 
be weakly reversible it is necessary that any such phantom reactant complex 
at least appear in the set of product reactants of M . 

4 Conclusion 

In this paper we have summarized what are, to the best of our knowledge, 
the most comprehensive results on the topic of chemical reaction networks 
with qualitatively equivalent dynamics under the assumption of mass-action 
kinetics. We have attempted to unify these results into a common theory 
which we have called conjugate chemical reaction network theory. We have 



of N is 



6t 3 = kik 2 k 3 . 
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presented a result (Theorem I3.2p which goes beyond the results currently 
in the literature and shown in Section 13.31 how this result can be used to 
expand the scope of weakly reversible networks theory. 

There are several immediate avenues for future work in the field of con- 
jugate chemical reaction networks. 

1. It is an open question the extent to which non-linear mappings h : 
I£>o i — y RPq can offer insight into the behaviour of chemical reaction 
networks. In order to preserve mass- action dynamics, it is clear that 
we should select h(x) from the class of polynomial transformations. 
The algebra, however, is non-trivial except for simple cases, and it is 
known that non-linear mappings need not preserve the form of the 
reactant complexes or all aspects of the qualitative behaviour of the 
mass-action system (see Appendix). 

2. The primary focus of Section 13.31 was determining when a reaction 
network N could be shown to be conjugate to a weakly reversible 
network A/ 7 . For all but the smallest examples, however, it is typically 
not feasible to check all of the conditions required for conjugacy over 
the set of possible weakly reversible target networks by hand. The 
development of algorithms and computer software which can efficiently 
check for viable weakly reversible target networks A/ 7 , and which makes 
full use of the techniques of "splitting" reactions and using "phantom" 
reactions, is a primary interest. 

We note that Szederkenyi and Hangos have presented some computer 
algorithms for the related problem of determining when a network 
has mass-action kinetics (|3j) equivalent to that of a complex balanced 
system for a specified set of rate constants |30| . 

3. We have assumed throughout this paper that a chemical reaction net- 
work J\f obeys the mass-action kinetics ([3]) , which corresponds to 



with Ri(x) = feillfci^/' Within the literature, and in particular 
within the systems biology literature, one often encounters other ki- 
netics. A popular choice for an enzyme-catalysed single-substrate re- 
action is the Michaelis-Menten kinetics 



dx 
~dt 



r 



^-zOi^(x) 



K M + x 
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where x is the concentration of the substrate. Another popular choice 
is a Hill function [T7]. Some authors consider general properties on the 
i?i(x), such as convexity, strict monotonicity with respect to reactant 
concentrations, etc. [3l[25]. It remains an open question as to how 
conjugacy can be applied to these alternate kinetic schemes and which 
aspects of qualitative dynamics are preserved. 

Appendix 

Comparing the qualitative behaviour of two conjugate mass-action sys- 
tems is more challenging when the diffeomorphism h : M™ i— >■ R™ is non- 
linear. In general, many qualitative properties are not preserved by non- 
linear mappings. 

For example, consider the mass-action networks Af and Af' given by 

Af : A O 



and 

Af' : A A 2A. 

Solutions of ([3]) associated with Af are given by x(t) = xoe~ kt while solutions 
of ([3]) for Af' are given by y(t) = yoe kt which can be mapped into one 
another via the non-linear mapping y(t) = l/x(t). It is clear the qualitative 
behaviour of Af is not preserved through the conjugacy mapping. 

More complicated disparities in behaviour can occur in higher dimen- 
sions. For instance, the reaction network Af given by 

Af: Ai A A 2 

is conjugate to the reaction network Af' given by 



Af' 



via the mapping 



Ai A O 

A 2 A 2A 2 

2Ai +A 2 A Ai + A 2 

A x + 2A 2 A Ai + 3A 2 



yi = x\ x~ 2 
V2 = ^r 1 3 
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In this case, we have the limiting behaviour x\(t) — > and x^ft) — > x\(0) + 
^2(0) for M but yi(t) — > and y2(t) — > 00 for M' . The qualitative behaviour 
is significantly different. 

There are still examples, however, for which non-linear conjugacies are 
capable of demonstrating qualitative equivalence of solutions. For example, 
consider the network 

M : A*± 2A. 

It is easy to see that M is conjugate to the reaction networks Mb given by 

ki/b 

M b : A +± (1 + b)A 

k 2 /b 

for b G Z>o via the diffeomorphism y{t) = xit) 1 ^ where y(t) is the solution 
of the mass-action system §3§ for Mb- The set of reactant complexes changes 
in all cases except 6 = 1; however, the qualitative behaviour of the systems 
are identical. 
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